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ABSTRACT 

We use A^-body simulations of star clusters to investigate the possible dynamical ori- 
gins of the observed spread in core radius among intermediate-age and old star clusters 
in the Large Magellanic Cloud (LMC). Two effects are considered, a time- varying ex- 
ternal tidal field and variations in primordial hard binary fraction. Simulations of 
clusters orbiting a point-mass galaxy show similar core radius evolution for clusters 
on both circular and elliptical orbits and we therefore conclude that the tidal field 
of the LMC has not yet significantly influenced the evolution of the intermediate-age 
clusters. The presence of large numbers of hard primordial binaries in a cluster leads 
to core radius expansion; however, the magnitude of the effect is insufficient to explain 
the observations. Further, the range of binary fractions required to produce significant 
core radius growth is inconsistent with the observational evidence that all the LMC 
clusters have similar stellar luminosity functions. 
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1 INTRODUCTION 

Star clusters are vital laboratories for testi ng theories of stel- 
lar d ynamics and stellar evolution (e.g. iMevlan fc Heggie] 
Il997l) . The globular clusters orbiting the Milky Way con- 
stitute a sample of about 150 evolved stellar systems with 
ages > 10 10 years. These provide useful constraints on the 
long term evolution of clusters but information about the 
paths which led to these end states can only be inferred via 
comparisons with the end points of numerical cluster simu- 
lations. In contrast, the system of rich star clusters in the 
Large Magellanic Cloud (LMC) includes clusters of compa- 
rable mass to Galactic globular clusters but which span the 
entire age range from less than 5 x 10 6 years (e.g. R 136) to 
about 10 10 years (e.g. NGC 2019). This allows us to observe 
massive star clusters at various epochs in their evolution 
and makes it possible to identify evolutionary trends, under 
the assumption that the external processes affecting cluster 
evolution have not changed significantly over the lifetime of 
the LMC. With this aim in mind, HST project GO7307 (e.g. 
iBeaulieu et"al"lll999D obtained deep images of four pairs of 
LMC clusters. Within each selected pair the clusters are 
approximately coeval, have similar masses and metallicities 
and lie at comparable radii from the LMC centre but dif- 
fer significantly in other properties such as core radius. The 
pairs range in age from about 10 7 years up to 10 10 years. 
The youngest clusters (NGC 1805 and NGC 1818) are old 
enough to have expelled all their residual gas and so are pure 
JV-body systems. 

Fig. shows a plot of the observed core radii of star 



clusters in the LMC as a function of their ages, with the clus- 
ters in the GO7307 HST sample identified by name. There 
is a clear trend that the spread in core radius is an increas- 
ing fmr^tionof cluste r age. This trend was first identified 
by lElson et alJ (Il989l) from ground-based imag es and was 
recently con firmed bv | Ma ckcv & Gilmore (2003) using HST 
archive data. lMackev fc Gilmoref <l2003l calculated core radii 
for 53 LMC clusters in a uniform manner with the aim of ob- 
taining an unbiased data set. They discuss in some detail the 
uncertainties in the data and conclude that selection effects 
ar e not responsible for the observed trend. They suggest, 
as lElson et alJ lll989l) had earlier argued, that the increased 
spread in core radius is an evolutionary phenomenon: all 
clusters formed with relatively small core radii (~ 2 — 3 pc) 
and subsequently some clusters experienced core expansion 
while others did not. 

lElson et alJ lll989l) discuss the possibility that inter- 
cluster variations in the slope of the initial mass function 
(IMF) could have produced the observed spread in core ra- 
dius. They point out that mass loss from the increased num- 
bers of massive stars found in clusters with IMFs flatter than 
the standard Salpeter IMF could lead to core expansion of 
th e observed ma gnitude on a timescale of about 10 8 ' 5 years. 
As lElsonl lll99lf) notes, however, the IMF slope required to 
explain the largest core radii would lead to disruption of the 
clusters on a timescale of a few x 10 7 years. More directly, 
based on a detailed obse rvational analysis of th e clusters in 
the GO7307 HST sample. ldeGriis et all J2002bl) have shown 
that the slopes of their luminosity functions are similar to 
each other within the observational uncertainties, strongly 
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suggesting that their IMFs were very similar down to stellar 
masses of about 0.8 — I.OMq. The absence of inter-cluster 
population variations implies that stellar evolutionary mass 
loss cannot be invoked to explain the observations. 

As a star cluster evolves, encounters lead to a redistribu- 
tion of energy among the stars. Heavier stars and binaries 
acquire lower velocities, on average, and sink towards the 
cluster core. This process of mass segregation could enhance 
the heating effect of a population of primordial binaries be- 
cause the binaries would tend to fall to the centre where en- 
counters occur more frequen tly. Observations of young clus- 
ters (e.g. lHillenbrandlll997t) and numerical simulations of 
star formation (e.g. IBonnell et alJ l2001,l show that so-called 
primordial mass se gregation can a l so resu lt from the cluster 
formation process. Ide Griis et all (j2002bT) have shown that 
there is strong evidence for mass segregation in all the young 
and intermediate-age clusters in the GO7307 sample. Fur- 
ther, the amount of mass segregation appears to be similar 
in all clusters and cannot therefore account for differences 
in the cluster core radii. 

The aim of this paper is to investigate, by means of N- 
body simulations, whether dynamical processes could give 
rise to the observed core radius-age relation. In this con- 
text, the pair of clusters NGC 1868 and NGC 1831 are of 
particular interest. These clusters have very similar ages and 
metallicities and t heir ma sses differ by at most a factor of 
two (see Table 1 of Ide Griis et all2002bh . However, as Fig.[T] 
shows, the core radius of NGC 1831 is almost three times 
larger than that of NGC 1868. Given that the cluster ages 
are only around 500 Myr any process which has affected the 
core radii must be capable of producing significant evolution 
even on this relatively short time scale. In this paper, we con- 
sider the dynamical effects of (1) variations in the external 
tidal field experienced by the clusters and (2) differences in 
the initial numbers of dynamically active binary stars in the 
clusters. Given that the crossing times at the half-mass radii 
of these clusters are about 10 7 yrs and the times to orbit the 
LMC are about the same order of magnitude, it is reason- 
able to suppose that either of these processes could have 
begun to produce observable effects on cluster structure on 
the relevant time scale. 

During the past decade, numerical simulations of star 
clusters have become increasingly realistic owing to dramatic 
increases in computational speeds and significant imp rove- 
ments in programming algorithms (e.g. Aarseth 1999). The 
JV-body code use d in this paper incorporate s the stellar evo- 
lution models of iHurlev etafl fcOOCl l2002h . This makes it 
possible to analyse the output from the simulations in an 
identical manner to that used for observational data, thus 
allowing meaningful confrontation of models with observa- 
tions. This is particularly important when discussing cluster 
properties such as core radius because the relation between 
the theoretical and observational definitions of this quantity 
depend on the conversion between the mass function and 
luminosity function of the cluster stars. For a real cluster 
this conversion is non-trivial because it depends on, among 
other quantities, the cluster metallicity and is further com- 
plic ated by the possible presence of mass segregation (see 
e.g. Ide Griis et al. 2002a). In this paper, all comparisons be- 
tween models and observations are made in terms of directly 
observable properties. 

In the interests of clarity, we include here a note on the 
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Figure 1. Observed core radii R c of LMC star clusters as a func- 
tion of age. The clusters in the GO7307 HST sample are iden- 
tified by name, and matched pairs are join ed by dotted lines. 
(Following Fig. 1 of Ide Griis et alj J2002rj) . data taken from 
iMackev fc Gilmord bOOfl l. 

definition of the core radius. As lCasertano fc HulJl|l985F) dis- 
cuss, theoreticians, numerical simulators and observers as- 
sign different meanings to the term "core radius" , reflecting 
the quantity of greatest relevance to their purpose. In theo- 
retical work, the core radius r c refers to the natural length 
scale of the system un der study, f or example the scale ra- 
dius of the widely used iKine! l)l96fil) profile. The operational 
definition of the core radius commonly used in discussions 
of iV-body simulations is that of the density radius ra- This 
is the density-weighted average of the distance of each star 
from the density centre (lAarsethl200llh where the local den- 
sity at each star is estimated from the mass within a sphere 
containing the star's six nearest neighbours and the density 
centre is defined as the density-weighted mean position of 
the stars. Observationally, since the available information is 
either in the form of cluster surface brightness profiles or lu- 
minosities of individual stars instead of mass densities and 
stellar masses, the core radius R c corresponds to the radius 
at which the projected surface brightness falls to half its 
central value. In the case of a cluster of equal-mass stars 
in equilibrium, there is generally a well-defined relat ion be- 
tween r c and R c - For example, for a lPlummeil (Il91lf> model 
it is straightforward to show that R c = 0.6435r c . It is im- 
portant to note, however, that during the evolution of the 
cluster this relation is not static and that there is no general 
relation between and either r c or R c . In this paper, we 
use the term "core radius" exclusively to refer to the ob- 
servational quantity R c . In Section l3,4l we will compare the 
evolution of and R c , in order to emphasize the care which 
must be taken in comparing observations with simulations. 

In Section |3 we describe briefly the iV-bodycode used 
in the simulations, as well as the pipeline used to calculate 
observational properties of the simulated clusters. Section [3] 



© 2003 RAS, MNRAS OOO.ITinU 



Core Radii of star clusters 3 



Table 1. Summary of JV-bodyruns performed. No, Mo and /j> q 
are the initial number of stars, total stellar mass and primordial 
fraction of hard binaries, respectively. The clusters were all sim- 
ulated in the tidal field of a point mass representing the LMC, 
and each simulation ran until time T max . 



Name 


No 


M (M s ) 


/b,0 


Orbit 


Tmax 


Circ 


5000 


2390 





Circular 


1.0 Gyr 


Ecc 


5000 


2410 





Elliptical 


1.0 Gyr 


Binary 


4500 


2310 


0.5 


Circular 


1.7 Gyr 


Binary2 


4500 


2390 


0.25 


Circular 


1.6 Gyr 


Binary3 


4500 


2180 


0.1 


Circular 


1.0 Gyr 


Binary Ecc 


4500 


2340 


0.5 


Elliptical 


1.6 Gyr 



considers the effects of time variation of an external tidal 
field on the evolution of the cluster core radius while Sec- 
tion l4.1l deals with clusters moving in a steady tidal field but 
with different initial binary fractions. In Section [4.21 we in- 
vestigate the combined effect of both mechanisms by study- 
ing binary-rich clusters in a time- varying tidal field and the 
results of these simulations are discussed in the context of 
the LMC clusters in Section 14.31 Finally in Section |S] we 
summarise our conclusions. 



2 A-BODY SIMULATIONS 

iV-body simulations are a very valuable tool in the study 
of star clusters because they allow many different pro- 
cesses to be studied with a minimum of simplifying as- 
sumptions. The simulations discus sed below were p erformed 
using the NBODY4 code (e.g. lAarsethl Il999f) running 
on the HARP-3 and G RAPE-6 special purpose comput- 
ers iMakino et alll997tl at the Institute of Astronomy, Cam- 
bridge and on a GRAPE-6 board at the American Mu- 
seum of Natural History. The NBODY4 code uses the 4th 
order Hermite integration scheme to follow the orbits of 
the clust er stars and utilises state-of-the- art regularisation 
schemes (iMikkola fc Aarsethlll993l . Il998l) for the accurate 
and efficient treatment of the internal motion of binary and 
higher-order subsystems. It also incorporates a full set of 
stellar evolution models, in the form of parameterized func- 
tions, which are used to follow t he evolution of single stars 
and binaries jHurlev et al.ll200lTl . Of particular importance 
to the present study is the treatment of mass loss. Any gas 
which is lost from an evolving star and is not accreted on to 
a binary companion is assumed to leave the cluster instanta- 
neously, thereby reducing the depth of the cluster potential 
well. A realistic treatment of mass loss is therefore essential 
to allow the study of the evolution of such cluster properties 
as the core radius. 

Within the A-bodycode the equations of motion are 
integrated in scaled units such that G — 1 and the to- 
tal m ass and virial radius of the cluster are also set to 
unity jHeggie fc Mathiedll986lh For a star cluster in virial 
equilibrium the initial energy in these units is —1/4 and 
the crossing time is 2\/2. Given the total mass of the clus- 
ter in solar masses and an appropriately chosen length scale 
(which determines the conversion from Af-body length units 
to physical units) it is straightforward to obtain the conver- 



sion factors for time and velocity from A?-body units to Myr 
and kms -1 , respectively. 

In order to facilitate direct comparison of the results of 
our simulations with the observed data, we determine the 
core radius for each model c luster using the same proc ess- 
ing pipeline as that used bv lMackev fc Gilmord i2003T) . In 
summary, the stellar data from the model clusters (in the 
form of position and luminosity information for each mem- 
ber star) are binned radially to obtain a proje cted surface 
brigh tness profile. The three-parameter model jElson et alJ 
1987) 

tir)=IM>(l + £) 2 (1) 

is then fitted to the cluster profile. Here, no is the central 
surface brightness, 7 determines the fall-off of the surface 
brightness profile at large radii and R c is related to the 
length scale a via Rc = a(2 2h - 1) 1/2 . Full details of the 
fitting procedure used an d discussion of the associa ted un- 
certainties are presented in lMackev fc Gilmord i2005t) . In or- 
der to reduce the numerical noise due to the small number of 
stars in our simulated clusters we observe each cluster along 
three perpendicular lines of sight and obtain the median 
core radius from these three estimates. This helps to reduce 
the impact of individual bright stars which might otherwise 
skew the core radius estimation. Further, each A-bodyrun 
is performed four times using different random realisations 
of the initial cluster model and the results are averaged in 
order to ensure the robustness of any observed trends. 

Initial ly, all the model clusters are generated according 
to a lPlummerl (| 19 1 lin density law 

P(^7^L (2) 

following the prescription of lAarseth et al 1 (119741) . Here rp is 
the Plummer length scale and po is the central mass density. 
Observations of the young clusters in t he LMC show that 
they do not appear to be tidally limited iMackev fc Gilmorel 
2003). In this case, the Plummer profile is a reasonable 
model to use for the initial density profile given that the 
fall-off in its outer regions lies between that of a tidally 
truncated King model and the observed profiles of the young 
LMC clusters. 

The i nitial distribution of stellar masses is drawn from 
the IMF of lKroupa et alJ (|l993h which is a three-part power 
law fit to stellar data in the solar neighbourhood. We take 
the lower and upper mass limits to be O.IMq and 5OM0, 
respectively. We assume solar metallicity for the stars in our 
clusters. This is an appropriate choice for the younger clus- 
ters - for example, Ijohnson et alJ (1200 ill have shown that 
NGC 1805 and NGC 1818, whose ages are about 10 and 
25 Myr respectively, are well matched by isochrones of ap- 
proximately solar metallicity. Since the pairs of clusters in 
the GO7307 HST sample were chosen to have similar metal- 
licity, we know a priori that the observed range of clus- 
ter core radii cannot be due to metallicity variations. We 
have therefore chosen to perform all our simulations at solar 
metallicity, and postpone discussion of the effects of metal- 
licity variations on star cluster evolution to future work. 

Table ^summarises the simulations which we have per- 
formed. The clusters in these runs initially contain 4500 or 
5000 stars and have a mean mass of about 2300 Mq. Clearly, 
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the LMC clusters are more massive and contain greater num- 
bers of stars than our model clusters and therefore direct 
comparison of the results of our simulations with the ob- 
served data must be approached with caution. However, the 
effects of increasing N in our simulations can be estimated 
by considering the various time scales in the problems we 
study. Detailed discussion of these scaling issues is presented 
in Sections 13.41 and 14.31 where we argue that none of our 
conclusions is strongly dependent on the small size of our 
simulated clusters. 



3 TIME- VARYING TIDAL FIELDS 
3.1 Model Initialisation 

We first consider the effects of the external tidal field 
of the LMC on the evolution of the cluster core radii. 
Theoretical arguments and numerical simulations strongly 
suggest that tides have played a significant role in the 
evolution of star clusters and c luster systems in galax- 
ies such as the Milky Way (e.g. [ Gnedin fc Ostrike:" 
Vesperini fe Heggielll997t IVesperiri1^)9STBaTmigaref' 



Fall & Zhang 2001) . In addition, observations of tidal tails 
around at least t wenty Milky Way globula r clusters (e.g. 
iLeon et al.l l200Ct lOdenkirchen et all 1200 fll provide com- 



pelling evidence of the ongoing disruption of star clusters by 
the Milky Way. The gravitational field of the LMC (mass « 
1O 1O M0) is much weaker than that of a g alaxy like the Milky 
Way (mass « 10 12 Mq). However, as iMackev fc Gilmord 
(2003) show (their Fig. 15), there is a deficit of low-mass 
clusters older than 10 10 years. If we assume that the cluster 
population originally contained a spectrum of cluster masses 
similar to that of the present-day young clusters, their ab- 
sence at late times strongly suggests that tides have been at 
work on the cluster system in exactly the man ner expected 
from the models of iGnedin fc Ostrikerl (Il997f) . Given that 
the GO7307 HST sample of 8 clusters was selected on the 
basis that they lie at similar distances from the LMC cen- 
tre but display a wide range of core radii, the steady tidal 
field experienced by a cluster on a roughly circular orbit 
cannot be the origin of the observed trends. Further, given 
the intrinsic weakness of the LMC field at a radius of 4kpc, 
where our clusters lie, one might not expect the tidal field to 
play a significant role in cluster evolution. We therefore con- 
sider clusters on non-circular orbits. Such clusters experience 
tidal shocks during their orbits and we investigate whether 
these impart sufficient energy to the cluster t o produce dif- 
ferences in the evolution of R c . lElsonl (12000^1 discusses the 
possible role of disk shocking in the internal evolution of the 
intermediate-age LMC clusters, but questions whether any 
of the clusters is on an appropriate disk-crossing orbit for 
this to be an important effect. In this section we consider 
instead the tidal shocking produced by a point mass galaxy. 

In order to investigate the role of the LMC tidal field in 
the evolution of its cluster system, it is necessary to incor- 
porate an appropriate external gravitational field into the 
simulations. We model the LMC potential as that of a point 
mass of 9xl0 9 M Q . This is a rather drastic oversimplification 
of the actual situation. We note, however, that this model 
exaggerates the effect of the tidal field, as any extended dis- 
tribution of mass will have a more slowly varying gravita- 
tional potential. Further, in the radius range 2 — 8 kpc, in 



which our model clusters orbit, the gradient of the potential 
is within a factor of two of that of t he detailed LMC mass 
model of Ivan der Marel et alJ J2002t> . Thus our simulations 
constitute a plausible upper limit to the possible effects of 
tides. 

In the iV-body simulations, the equations of motion 
are integrated in an accelerating but non-rotating reference 
frame centred on the centre of mass of the cluster and whose 
axes are aligned with fixed spatial directions. This choice of 
frame simplifies the equations of motion for a cluster on a 
non-circular orbit despite having less symmetry than the 
more usual rotating frame used for steady tidal field s whose 
axes always point towards the attracting centre llHeggiel 
l200lT) . Implementation of the time varying tidal field re- 
quires some care because mass loss from evolving stars or 
the escape of stars from the cluster causes the centre of 
mass to acquire a non-zero velocity relative to the frame of 
integration. It is therefore necessary to re-centre the frame 
of reference to coincide with the rest frame of the centre 
of mass following each such event. All internal and external 
forces must then be reinitialized to ensure consistency with 
the re-centred frame during subsequent integrations. 

We have performed simulations for clusters on elliptical 
orbits with perigalactic radii R p « 2 kpc and apogalactic 
radii J? a ~ 8 kpc, probing the region in which most of our 
sample clusters lie but largely avoiding the complications of 
the inne r regions of the LMC which may contain a stellar 
bar fe.g. lZhao fc Eva ns 2000). We also modelled clusters on 
circular orbits at a radius of 6 kpc which corresponds to the 
time-averaged radius of the elliptical orbits. This allows us 
to isolate the effect of the time variation of the tidal field 
because both sets of clusters experience the same mean tidal 
field. 

We calculate the initial perigalactic tidal radius r t , p in 
the approximation that the cluster is on a circular orbit at 
a radius R p . In this case, r t , p is determined from the cluster 
mass M ; the mass of the point mass galaxy Mq and R p 
via (e.g. lKindlT962T) 



rt, P 



Rr, 



M 



3GM G 



1/3 



(3) 



The apocentric tidal radius rt,a is determined in a similar 
way from J? a . The unit of length, which determines the con- 
version between iV-body length units and physical units, is 
then chosen such that the limiting radius of all clusters is 
initially equal to r t , p . However, at the start of the simula- 
tions the clusters are placed at the apocentre of their orbits 
about the LMC which means that they do not fill their tidal 
radii. If, instead, the clusters initially extended to r t , a and 
were not truncated at rt lP then the removal of stars outside 
rt, p during the first few perigalactic passages could mask 
any heating effects at radii interior to rt, p due to the time 
variation of the tidal field. In contrast to a King model, 
the Plummer density law © formally extends to infinity 
and so the outer radius of a cluster is not well-defined. The 
length unit we use yields a pericentric tidal radius close to 
the (random) position of the outermost star in our initial 
clusters. The internal stellar velocities are scaled to ensure 
virial equilibrium at the start of the simulation. 

The tidal radius of a cluster on an elliptical orbit does 
not have a unique value. It varies continuously both as the 
cluster moves along its orbit and when the cluster loses mass 
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through stellar evolution or the escape of member stars. 
However, for clusters whose internal relaxation times are 
longer than their orbital perio ds, tidal strippi ng tends to 
truncate the star cluster at r t , p iOh fc Lmlll992l) . These two 
time scales are of similar magnitude for our simulated clus- 
ters so we expect two-body relaxation processes to lead to 
some expansion of the cluster beyond rt, p . We take account 
of this by defining escaping stars to be those stars outside 
twice the apogalactic tidal radius r t , a - Stars whose orbits 
carry th em outside th is radius are generally lost to the clus- 
ter (e.g. iHeggidboOlf) . Stars are thus free to populate the 
region between r tjP and 2rt, a - 

The choice of length unit determines the physical den- 
sity of the cluster and hence the physical rate at which two- 
body relaxation processes affect the evolution of the cluster. 
Our aim is to investigate the influence of a time-varying 
external tidal field so it is important to ensure that all in- 
ternal evolutionary processes proceed in an identical manner 
for each of our simulated clusters. We therefore use the same 
length unit for the clusters on circular orbits as we chose for 
those on elliptical orbits. This means that the clusters on 
circular orbits also do not initially fill their tidal radii and 
so necessarily allows expansion of these clusters. However, 
since the circular orbits are placed at the time-averaged ra- 
dius of the elliptical orbits this expansion should on average 
be the same for both sets of clusters. Any differences in the 
rates of expansion are due to the time variation in the tidal 
field. 

The effect of steady tidal fields on star cluster evolution 
has been the subject of papers by a number of authors (e.g. 
[ Giersz fc Heggidll997t lAarseth fc Heggid 11995 IBaumgardtl 
l200ll : IPortegies Zwart et al]|2002T) . However, until recently 
the problem of clusters moving in time-varying tidal fields 
had rec eived relatively litt le attention from TV-body simu- 
lations. IBaumgardtl il998f) uses the results of some simula- 
tions of single-mass star clusters orbiting a point mass to 
determine the appropriate choice of escape radius for sim- 
ulatio ns of mass loss. More recently. iBaumgardt fc Makinol 
( 2003) presented an extensive survey of the dynamical evo- 
lution of star clusters in external tidal fields. Our investiga- 
tion is complementary to theirs because we are looking at 
the evol ution of the observable global structure of the clus- 
ter while IBaumgardt fc Makind i200^l focus their attention 
mostly on the evolution of the stellar mass function within 
the cluster. The results of all our simulations are presented 
in terms of quantities which can be compared directly with 
those calculated from the observed data, thereby enabling 
us to discuss how different dynamical processes impact on 
observations. 

We follow the evolution of the clusters for about 5i r h 
(where t. r h is the half-mass relaxation time) This corresponds 
to Tmax ~ 1 Gyr in most cases. By this time, the clusters on 
elliptical orbits have lost more than 50 per cent of their stars 
making estimates of properties such as the core radius in- 
creasingly noisy. At regular intervals during the evolution, 
the data from the clusters are recorded and all relevant prop- 
erties of the clusters are calculated. 



3.2 Results 

The solid curve in Fig. [5] presents the evolution of the core 
radius for clusters on circular orbits about the LMC at a 
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Figure 2. Core radius evolution for iV-body simulations of clus- 
ters on orbits about a point mass representing the LMC. See text 
for discussion. 



distance of about 6kpc from the centre. The broken curve 
shows the evolution of R c for clusters on elliptical orbits 
with the same time-averaged galactocentric radius. It is ap- 
parent from this figure that the core radius evolution in both 
cases is indistinguishable, especially in view of the noise as- 
sociated with the output of small- TV simulations. During the 
first 10 8 yrs, the core radius is seen to decrease slightly from 
its initial value. This is because the calculation of Rc is dom- 
inated by the light of the most massive stars which rapidly 
become more centrally concentrated due to the effects of 
mass segregation. As Fig. |5| shows, after 10.5 Myr the half- 
mass radius of the more massive stars is already 30 per cent 
smaller than that of the less massive stars. 

After its initial contraction phase, the core radius be- 
gins to increase. This occurs because the clusters do not fill 
the tidal radius corresponding to their galactocentric posi- 
tions at the start of the simulations. A combination of mass 
loss from the rapid evolution of massive stars and two-body 
relaxation can therefore drive the overall expansion of the 
cluster seen in Fig. [3] The top panel of Fig. 2] shows that 
the cluster's half-mass radius r\ is increasing with time (see 
fourth panel of Fig. for long-term behaviour). During 
this period its rd (not shown) is roughly constant at about 
0.3 pc. 

The expansion of is driven by a combination of two 
processes which are illustrated in the lower panels of Fig. [I] 
The middle panel shows the evolution of the virial ratio Q. 
Initially the cluster is in equilibrium with Q = 0.5. After 
about 10 Myr, the cluster is temporarily displaced from equi- 
librium (Q > 0.5) by mass loss from the most massive stars. 
This is accompanied by the increase in the half-mass radius 
of the cluster shown in the top panel of Fig. [I] It is interest- 
ing to note that the Rc of the cluster does not change as a 
result of this mass loss. Stellar evolution cannot account for 
the sustained expansion of the cluster, however. The cluster 
loses about 5 per cent of its mass as a result of the above 
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Figure 3. Cumulative mass profiles for stars above (solid curves) 
and below (broken curves) the median stellar mass for a simulated 
cluster containing only single stars and on a circular orbit about 
the LMC. Profiles at two different times are shown: T = (upper 
curves) and T = 10.5 Myr (lower curves). There is clear evidence 
for mass segregation developing in the cluster. 



mass loss events. Using the virial theorem it is straightfor- 
ward to show that the gravitational radius R' of a stellar 
system with initial radius R and mass M which loses an 
amount of mass m is given by ( Hilisl lT980l) 



R' 



If M — m 



2 \\M- 



R. 



(4) 



A 5 per cent mass loss should lead to an increase of less than 
10 per cent in the size of the cluster which is insufficient to 
account for the cluster expansion seen in Figs. H and H 

The bottom panel of Fig. [I] shows the evolu- 
tion of th e velocity anisotropy p arameter /9 = 1 — 
(vl) /2{Vr) llBinnev fc Tremainel 19871) inside the 10 per cent 
Lagrangian radius of the cluster (the radius containing 10 
per cent of the total mass of the cluster) and in the outer 
parts (between the 80 per cent and 90 per cent Lagrangian 
radii). Initially, the velocity distribution throughout the 
cluster is isotropic (/3 = 0). While the central parts of the 
cluster retain this isotropy, in the outer parts the distribu- 
tion becomes increasingly radially anisotropic (/3 > 0). This 
indicates that the continued growth of is being fueled 
by encounters in the inner parts of the cluster which pro- 
duce a population of stars on extended radial orbits. The 
velocity distr ibution in the outer par ts of isolated clusters 
simulated bv lBaumgardt et alJ (I2002T) exhibited similar be- 
haviour during post core-collapse evolution; those authors 
similarly concluded that the growth of the outer Lagrangian 
radii at late times was due to two-body processes in the 
inner regions. 

Given the similarity between the core radius evolution 
of clusters on circular and elliptical orbits, we now proceed 
to investigate the extent to which the clusters on elliptical 
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Figure 4. Early evolution of the velocity structure of a cluster on 
a circular orbit. Top: evolution of the half-mass radius (solid 
curve) and core radius R c (broken curve); Middle: plot of the 
virial ratio Q (whose equilibrium value is 0.5); Bottom: velocity 
anisotropy parameter j3 of stars within the inner 10 per cent La- 
grangian radius (broken curve) and between the 80 per cent and 
90 per cent radii (solid curve). See text for discussion. 



orbits have been affected by the variation in the tidal field. 
Fig. presents the evolution of certain key internal proper- 
ties of two typical clusters, one on a circular orbit and one 
on an elliptical orbit. The top panel shows the galactocentric 
radius J? g along the elliptical orbit to emphasize the clear 
correlation between perigalactic passages and changes in the 
internal properties of the cluster on the elliptical orbit. The 
second panel shows that, as expected, during a perigalactic 
passage the rms velocity u rms of the stars in the cluster in- 
creases sharply. This velocity kick ejects the most energetic 
stars from the cluster. As a result of this mass loss, v ims 
settles to a value slightly below that of the cluster on the 
circular orbit once it returns to virial equilibrium. 

The third and fourth panels of Fig. show how t r h and 
rh evolve during the cluster orbit. The most noticeable effect 
is that the relaxation time for the cluster on the elliptical 
orbit is significantly reduced after the second perigalactic 
passage, so that the cluster evolves more rapidly following 
this encounter. This is due both to the decrease in the to- 
tal number of stars in the cluster and to the increase in the 
mean stellar mass m (shown in the bottom panel of Fig.0. 
The increase in m occurs because, on average, low-mass stars 
acquire larger velocities during pericentre passages and are 
hence more likely to escape from the cluster. The value of 
rh is seen to increase during the perigalactic passage but 
rapidly settles back to a value slightly below the circular- 
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orbit value. Thus while the tidal shocks experienced by the 
clusters on elliptical orbits do not affect the evolution of R c , 
they are nonetheless strong enough to modify other impor- 
tant properties of the cluster. As Fig. [2] shows, the impact 
of the shocks on the cluster properties increases with time - 
this is due to the decreasing number of bound stars in the 
cluster. 



3.3 Strength of tidal shocks 

Before drawing our final conclusions based on these simu- 
lations, we now investigate whether the lack of impact on 
R c is due to the weakness of the shocks we have simulated 
or whether the value of R c is intrinsically robust to tidal 
shocking. The strength of a tidal shock is generally presented 
in terms of t he dimensionless parameter /3 sh given by (e.g. 
ISpitzerlfl987l) 



8 



tsh ' 



(5) 



where t s h is an estimate of the destruction time of the cluster 
owing to energy input from the gravitational shocks. The 
time scale t s h is given by 



* sh = dETJdi' (6) 

where Eh is the energy of the cluster evaluated at the half- 
mass radius and dEh/dt is the rate at which stars at the 
half-mass radius acquire energy due to the shock. In a simu- 
lation containing a finite number of stars, evaluation of this 
expression requires some care. Averaging the energy changes 
of stars near the half-mass radius does not produce robust 
results because individual stars may gain and lose energy 
during the perigalactic passage through other means than 
interaction with the external tidal field. Instead we calcu- 
late the total energy binding the cluster as a function of 
time and estimate the energy input to the cluster during 
perigalactic passages. Fig. [S] compares the time evolution of 
the binding energy Eh of a cluster on a circular orbit with 
that of a cluster on an elliptical orbit. For a cluster in an 
exter nal potential Eh is defined as (e.g. lFukushige fc Heggid 
1995) 



E b = T + V + E T 



(7) 



where T and V are the total kinetic and potential energy of 
the cluster, respectively, and Et is the energy in the tidal 
field. The binding energy must be negative for the cluster 
to remain bound. The effects of the external field are visible 
as dramatic changes in Eh at about 180 and 520 Myr. Note 
that the sudden, temporary change in Eh for the circular run 
at about 700 Myr is due to a collision between two binaries 
that results in the the break-up of one of the binaries and 
the subsequent ejection of its components. The ratio of the 
energy produced by this event to the binding energy of the 
cluster emphasises the relatively fragile nature of the clusters 
we are simulating. 

During perigalactic passages, some stars acquire veloc- 
ities above the escape velocity of the cluster and this leads 
to an enhanced mass- loss rate (see Fig. 0. However, the fi- 
nite time required for such stars to escape means that Eh can 
temporarily become positive because it is calculated by sum- 
mation over all stars within twice the tidal radius of the clus- 
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Figure 5. (Top panel) Variation of galactocentric radius i? g with 
time for a cluster on an elliptical orbit (broken curve). The solid 
line shows the position of the circular orbit of equal time-averaged 
radius. (Lower four panels) Evolution of rms velocity » rms , half- 
mass relaxation time t T h, half- mass radius and mean stellar 
mass m. Curves for two cases are shown: (1) a cluster on a circular 
orbit (solid curves) and (2) a cluster on an elliptical orbit (broken 
curves). Note that the half-mass radius is calculated using all stars 
within a radius of 2rt, a °f the cluster centre. In the simulations 
all such stars are classified as non-escapers, although at any given 
time a fraction of them are formally unbound. 



ter. After the cluster has passed pericentre, stars continue to 
leave across the tidal boundary and Eh returns to a negative 
value indicating that the remaining cluster is still gravita- 
tionally bound. Following each shock, the cluster undergoes 
a revirialisation process during which the energy acquired 
from the tidal interaction is redistributed among the cluster 



© 2003 RAS, MNRAS 000.1111131 



8 M. I. Wilkinson et al. 



-0.1 



-0.2 



-0.3 




100 



200 400 600 
Time (Myr) 



800 1 000 



Figure 6. Evolution of the energy binding the cluster E^. Curves 
for two cases are shown: (1) a cluster on a circular orbit (solid 
curve) and (2) a cluster on an elliptical orbit (broken curve). En- 
ergies are given in standard iV-body units (G = Mo = l,Eb = 
-1/4). 
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Figure 7. Evolution of total cluster mass M for clusters on cir- 
cular (solid curve) and elliptical (broken curve) orbits in units 
of the initial mass Mo. Enhanced mass loss during the period of 
revirialisation which follows each perigalactic passage is seen in 
the elliptical orbit case. 



stars. This process results in further loss of mass from the 
cluster because, on average, the revirialisation process gives 
velocity kicks to stars which are three ti mes larger than those 
arisin g from the initial shock (see e.g. iGnedin fc Ostrikerl 
1999). 
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Figure 8. Evolution of the Lagrangian radii. Curves for two cases 
are shown (1) for a cluster on a circular orbit (solid curves) and (2) 
for a cluster on an elliptical orbit (broken curves). From bottom 
to top, the curves show the radii containing 10, 30, 50, 70 and 90 
per cent of the cluster mass. The heating effect from the external 
tidal field is evident, even at the radius containing only 30 per 
cent of the cluster mass. 



We take the energy input by the shock to be the dif- 
ference between the energy of the cluster after this reviri- 
alisation process is complete and the energy of the cluster 
on the circular orbit at this time (about T = 540 Myr). 
The time between the shocks is 340 Myr, which leads to 
t s h ~ 6800 Myr. At this time i r h for the cluster is about 
170 Myr, so the shock strength is about 0.0 2. Comparing this 
value with Fig. 13 of IGnedin et al.l (1 19991) and noting that 
the initial concentrations c = log 10 (rt/r<i) of our clusters are 
about 1.1 (where we assume that r t is the initial outer limit 
of the clusters), we find that our clusters lie approximately 
on the borderline of the region where shocks are important. 
We conclude that the shocks our clusters are experiencing 
would significantly affect their evolution on long time scales. 
Indeed, as Fig.|H|illustrates, even the 30 per cent Lagrangian 
radius of the cluster is affected by the tidal shocks, indicat- 
ing energy input into the cluster over a wide range of radii. 
Nevertheless, on a timescale of up to 10 9 yrs, the effects of 
tidal heating are insufficient to modify the observed R c of 
the clusters significantly. 



3.4 Discussion 

The top panel of Fig.0]illustrates an important point which 
is worth emphasizing in the context of comparisons between 
simulated data and observations, namely that R c and 
are not, in general, related in a simple way. In particular 
circumstances, for example at the start of our simulations 
where the model clusters are in strict dynamical equilibrium, 
it is easy to show that, for a Plummer model composed of 
equal mass stars, the two quantities are related by r\ — 
2R C . However, Fig.|l]shows that this relation quickly breaks 
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Figure 9. Comparison of the evolution of the density radius rd 
(broken curve) with the observational core radius R c (solid curve) 
for a single TV-body cluster. Top: Un-smoothed results from the 
simulations. Middle: Curves smoothed by averaging over five time 
bins. Bottom: Ratio of R c to r d calculated using the smoothed 
curves from the middle panel. See text for discussion. 



down as the cluster evolves away from its initial conditions. 
Thus, in the absence of further information, observational 
estimation of R c must not be assumed to imply knowledge 
of the mass distribution in the cluster. 

The same is true of the density radius defined in Sec- 
tion0 Fig.[5]shows the evolution of R c (solid curve) and rd 
(broken curve) for one of the clusters on an elliptical orbit. 
As the top panel of the figure shows, in a small- TV simula- 
tion, both 7i c and r<j are intrinsically noisy quantities: this 
is emphasised by the middle panel where we have smoothed 
the curves from the top panel using a window function of 
width five time bins. The bottom panel of Fig. [5] which 
shows the ratio of R c to rd for the smoothed data, demon- 
strates that the relation between the two quantities is not a 
simple proportionality and that the evolutionary trends in 
the two quantities can be quite different. The complicated 
relation between light and mass in star clusters is the result 
of a variety of processes such as stellar evolution and mass 
segregation which modify the mass and light distributions 
in different ways. Given the absence of simple relations be- 
tween the different core radii, great care must be exercised 
when we interpret observational data through comparison 
with numerical simulations or theoretical models. 

The simulations we have performed contain relatively 
low numbers of stars and are smaller and of lower total mass 
than the actual LMC clusters. However, the mean densities 



of our clusters are very similar to those of the LMC clus- 
ters and it is this which determines the extent to which an 
external tidal field can influence the internal dynamics of a 
cluster. On this basis, we can be confident that conclusions 
based on our simulations are relevant to the clusters in the 
LMC. Our simulations have shown that, on a timescale of 
1 Gyr, which corresponds to several t r h for our model clus- 
ters, tidal forces are not able to produce a spread in core 
radius similar to that observed among the LMC clusters. 

There are two obvious avenues for building on the sim- 
ulations we have performed. First, simulations of clusters 
containing a larger numbers of stars N could be performed. 
This would allow us to achieve cluster masses comparable to 
those of the real LMC clusters. However, increasing N in a 
cluster while holding the physical density of the cluster con- 
stant leads to an increase in i r h- The expansion seen in R c 
for our simulated clusters is driven primarily by two-body 
processes (which proceed on a time scale t r h) rather than 
by tidal heating and it will therefore be slowed by increas- 
ing N. Further, the disruptive effect of pericentre passages 
depends on the ability of a cluster to re-fill the region out- 
side its pericentric tidal radius r t . p during the time between 
close passages. As i r h increases this re-filling is reduced. The 
cluster will therefore have fewer stars outside rt lP and hence 
will be l ess affected by shocks than small-TV clusters (see 
Fig. 2 of lBaumeardt fc Makinol fcOOSt ). A cluster on an el- 
liptical orbit and which, as in our simulations, initially filled 
its pericentric tidal radius would remain truncated at this 
radius throughout its evolution because there would not be 
sufficient time between perigalactic passages for the region 
outside r t , p to be re- filled l|Oh fc Lirj|l992T) . Thus shocks 
should have less impact on the 7? c of larger- N clusters and 
the observed expansions of clusters on both circular and el- 
liptical orbits should remain similar. 

The second possible improvement would be to use a 
more realistic model of the LMC poten tial, for example one 
based upon the recent mass models of Ivan der Marel et alJ 
(2002). However, the actual LMC mass distribution is more 
extended than the Keplerian potential we have been us- 
ing in our simulations and therefore the gradients in its 
gravitational field are smaller. Consequently, the associated 
strength of its tidal effects on clusters moving on elliptical 
orbits is reduced. 

Thus, we conclude that neither of the above modifica- 
tions to our simulations would tend to increase the effect of 
the tidal field. The fact that the time-varying tidal field in 
our simulations did not result in changes in the evolution of 
R c in the extreme case of a low-mass cluster in orbit about 
a point-mass galaxy leads us to conclude that tides have 
not yet had time to play a significant role in the internal 
evolution of the intermediate-age LMC clusters. 

4 MODELS WITH PRIMORDIAL BINARIES 
4.1 Steady tidal field case 

In this section we present results for models moving in a 
steady tidal field but with different numbers of primordial 
binaries, in order to isolate the impact of large binary frac- 
tions on the evolution of the cluster properties. The clusters 
discussed here move on the same orbits as the circular orbit 
runs of Section [!J] 
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Based on an analysis of the observed bin a ry seq uence 
in the young cluster NGC 1818, lElson et all <1998h have 
shown that the binary fraction /b in the centre of this clus- 
ter is about 0.35. Here / b = N h / (N s + N b ) , where N s and jV b 
are the numbers of single stars and binary systems, respec- 
tively. To date, the binary fractions of the older clusters in 
our sample have not been determined. Thus, at present the 
binary fractio ns of t hese clusters are essentially free parame- 
ters. Heggie's 1^975) law implies that only hard binaries are 
dynamically important in the evolution of a star cluster. A 
binary is said to be hard when the magnitude of its binding 
energ y |-E b j n | exce eds the kinetic energy of a typical cluster 
star llHeggiall975l) . In the case of a binary with component 
masses mi and 7712 at a separation a moving in a cluster 
with mean mass m and velocity dispersion a the hardness 
condition is 



Gm\m2 1 2 

> —mo . 

2a 2 



(8) 



(We note that the NBODY4 code uses a dispersion of y/2a in 
this equation to determine the hardness limit for binaries, in 
order to take account of the larger velocity dispersion in the 
core of a centrally concentrated cluster.) In our simulations 
the boundary between hard and soft binaries occurs at a sep- 
aration of about 120 AU for binaries containing stars of the 
mean stellar mass in the cluster. Hard binaries act as heat 
sources because they become more tightly bound as a result 
of interactions with single stars a nd thereby impart kinetic 
energy to the stars they encounter. iHeggie fc Aarsethl dl992l) 
found that the density radii of clusters containing primor- 
dial binaries experienced shallower contraction during core 
collapse and subsequently tended to be larger than those of 
clusters containing only single stars. This result, based on 
simulations of single-mass clusters with binary fractions of 
about 0.03, motivates our investigation of clusters contain- 
ing large numbers of primordial binaries. Hereinafter, we will 
use /b to refer to the fraction of hard binaries in the cluster. 

As an extreme case, we consider first clusters with initial 
/ b = 0.5 and with the hard binaries randomly distributed 
throughout the cluster. In each binar y the mass of the pri - 
mary star is drawn from the IMF of iKroupa et alj <ll993t) 
and the mass of its companion is obtained by randomly 
choosing a mass ratio for the binary from a uniform dis- 
tribution constrained by the single star mass lim its (for a 
more detailed discussion, see iHurlev et al ] l200lft . The bi- 
n ary orbital parameters are generated from the distribution 
of lEggleton et aD (ll98ST) . The clusters are otherwise gener- 
ated exactly as in the single star runs of Section |3] We note 
that the inclusion of a large number of hard primordial bi- 
naries dramatically increases the computation time required 
for each simulation. 

Fig. 1101 compares the evolution of R c for clusters con- 
taining N = 5000 single stars and no primordial binaries 
with that for clusters containing a significant initial binary 
population. An increased binary fraction produces a notice- 
able effect on the evolution of R c . After 1 Gyr, the core radii 
of clusters containing 50 per cent primordial binaries are 
about a factor of two larger than those of clusters with no 
binaries. 

The results presented in Fig. 1111 give some insight into 
the mechanism by which a large binary fraction can drive 
an expansion of R c - The top panel shows that the rate at 
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Figure 10. Core radius evolution for TV-body clusters containing 
single stars only (short-dashed curve) and clusters containing 10 
per cent, 25 per cent and 50 per cent hard primordial binaries 
(short- long dashed, long-dashed and solid curves, respectively). 
All clusters are on circular orbits about the LMC. See text for 
discussion. 



which mass is lost from a cluster with initial /b = 0.5 is 
significantly greater than in the case of a cluster with no 
initial binary population. After 1 Gyr this cluster has lost 
almost 25 per cent more mass than the cluster containing 
only single stars. The evolution of the energy binding the 
cluster E b is shown in the second panel of Fig. 1111 In a clus- 
ter with a large initial binary fraction, 73b rapidly attains 
a value significantly higher than that of a cluster of single 
stars. The enhanced mass-loss rate due to the presence of 
the binaries decreases the depth of the cluster's potential- 
well; the magnitude of the binding energy is further reduced 
by the heating of the remaining bound stars through weak 
encounters. As Fig. 1121 shows, the energy input due to bi- 
nary heating is distributed throughout the cluster. All the 
Lagrangian radii of the cluster rapidly attain larger values 
when the cluster contains primordial binaries. We note that, 
as was discussed in Section 13.21 the growth of the clusters 
containing only single stars is due to the fact that the clus- 
ters do not fill their tidal radii at the start of the simulations. 
We use identical initial conditions for the clusters contain- 
ing binaries - the excess heating observed in those clusters 
is due to the presence of the hard binary population. 

The third panel of Fig. 1111 shows the internal binding 
energy of the primordial hard binaries -Ebm as a function of 
time. In calculating £" b in we include the binding energy of 
those hard binaries which have escaped from the cluster be- 
cause we are interested in the rate at which encounters lead 
to the hardening of binaries. This plot demonstrates that, 
as expected, the hard binaries become harder with time and 
so impart energy to the cluster. It is well known that en- 
ergy generation by hard binaries halts core collapse (e.g. 
IHeggie fc Aarsethl Il992'l . However, Fig. 1111 demonstrates 
that even in the very early stages of a cluster's evolution, the 
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Figure 11. Comparison of the evolution of the total cluster mass 
M (top panel) and the energy binding the cluster Eb (middle 
panel). Two cases are shown (1) a cluster without primordial 
binaries (broken curves) and (2) a cluster containing 50 per cent 
primordial binaries (solid curves). The total mass is normalised 
by its initial value Mo and Eb is given in standard TV-body units. 
The bottom panel shows the evolution of the internal binding 
energy of the primordial binaries Eb\ n (in standard TV-body units) 
in the run containing primordial binaries. The decrease in -Ebin 
is the result of interaction-induced binary hardening. Note that 
Ebin includes the internal energy of escaped binaries. See text for 
discussion. 



heating effect of a large primordial hard binary population 
can have a significant impact on the energy budget of the 
cluster. 

We note that the occasional large changes in £"bin are 
the result of stellar evolutionary processes and therefore do 
not contribute to the direct heating of the cluster. For ex- 
ample, at T — 34 Myr the components of one binary system 
merge to form a single star: prior to the merger the binary is 
hardened by tidal interactions between its components. The 
massive star thus produced explodes as a supernova soon af- 
ter and the velocity kick received by the remnant as a result 
of the explosion is reflected in the spike in Eb- The promi- 
nence of this event is due to the relatively small number of 
stars in the simulation. 

Fig. 1101 shows that the R c evolution of clusters with 
fb = 0.25 is very similar to that of clusters with fb = 0.5 
initially. This suggests that the heating effect of primordial 
binaries saturates above a certain value of the binary frac- 
tion /b.crit- For larger values of fb, the R c evolution is un- 
changed, although the mass-loss rate is greater for increased 
fb- From Fig. HOI we see that /b,crit for the clusters we have 
simulated appears to lie in the range 0.1 — 0.25 because an 
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Figure 12. Evolution of the Lagrangian radii for clusters with 50 
per cent primordial binaries (solid curves) and with single stars 
only (broken curves). From bottom to top, the curves show the 
radii containing 10, 30, 50, 70 and 90 per cent of the cluster mass. 
The heating effect of the primordial binaries is visible at all radii. 



initial binary fraction of 10 per cent has a noticeably smaller 
impact on R c than one of 25 per cent. We are currently inves- 
tigating this effect in more detail, including its dependence 
on the number of stars in the simulation: the results will be 
presented elsewhere. Here we note only that the saturation 
may arise from the fact that at early times only a relatively 
small number of binaries take part in the heating process. 
In our simulations only about three per cent of the binaries 
are dynamically hardened during the first 100 Myr. We find 
that this fraction appears to be independent of the value of 
fb- Whether a particular binary experiences an encounter 
which leads to hardening depends on the masses and sepa- 
ration of the components as well as the local stellar density 
in the neighbourhood of the binary. For a fixed binary mass, 
however, the hardest binaries have the smallest separations 
and hence are least likely to interact with other stars. This is 
reflected in our simulations: for a given binary mass, all the 
binaries which undergo early hardening are initially among 
the more weakly bound pairs. 



4.2 Time-varying tidal field case 

In our final set of simulations we examine clusters contain- 
ing large numbers of primordial binaries and moving on non- 
circular orbits about the LMC. We have considered the ex- 
treme case of fb = 0.5 but, as was shown in the previous 
section, these results also apply to somewhat smaller values 
of fb- As Fig. 1131 shows, the core radius of a cluster on an 
elliptical orbit and containing primordial binaries exhibits 
stronger growth than that of a cluster of single stars on a 
circular orbit: after about 600 Myr the core radius of the 
former cluster is almost a factor of two larger. The increase 
in core radius is solely due to the presence of binaries. As in 
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Figure 13. Core radius evolution for JV-body simulations of clus- 
ters on different orbits about a point mass representing the LMC 
with (solid and long-dashed curves) and without (short-dashed 
curve) hard primordial binaries. The solid curve represents a clus- 
ter on an elliptical orbit, while both broken curves represent clus- 
ters on circular orbits. See text for discussion. 

the case of clusters containing only single stars (Section |3J, 
a cluster with the same primordial binary fraction but on a 
circular orbit exhibits almost identical core radius evolution. 
However, the increased mass- loss rate which results from bi- 
nary heating leaves the cluster on the elliptical orbit more 
vulnerable to tidal disruption. After only 600 Myr, the clus- 
ter has lost almost 75 per cent of its mass and the core radius 
determination is consequently very uncertain. By 1 Gyr it is 
virtually disrupted and retains only about 10 per cent of 
its initial mass while the circular orbit cluster still contains 
about 45 per cent of its initial mass. 

4.3 Discussion 

Our simulations demonstrate that the presence of a large 
population of hard binaries can produce significant ex- 
pansion of the core radius of a cluster. In agreement 
with many previous studies (e. g. |SDitze^^^Matjiieul ll98Ct 
iHeggie fc Aarsethll99 Aide La Fuente Marcosil996l) we have 
thus shown that binaries constitute a dynamically signifi- 
cant population in star clusters. In particular, we have found 
that primordial binaries can rapidly affect observable cluster 
properties such as R c . In the context of the LMC clusters, 
however, it is clear that the presence of binaries alone cannot 
account for the the observed trends in R c and in particular 
for the factor of three difference in core radius between NGC 
1831 and NGC 1868. Even in the extreme case in which one 
cluster has no primordial binaries at all while the other clus- 
ter has a hard binary fraction of 25 — 50 per cent, the differ- 
ence in R c over the lifetime of the intermediate-age clusters 
is at most a factor of two. The core radius of a cluster con- 
taining 10 per cent binaries expands by only a factor of 1.5 
over the same time scale. Further, the similarity of the ob- 



served present-day lu minosity functions of the clusters in 
the GO7307 sample dde Griis et al.ll2002bl) argues against 
large inter-cluster variations in the IMF and hence makes 
significant differences in binary fraction unlikely. 

The initial radial distribution of hard binaries has an 
impact on their heating effect on the cluster stars (e.g. 
Ide La Fuente MarcodlT996l) and may also have a bearing on 
the saturation of the heating effect which we observed in 
Section [4. II A binary population which was more centrally 
concentrated than the single stars as a result of either pri- 
mordial or dynamical mass segregation would be expected 
to produce more heating because of the increased interac- 
tion rates in the inner parts of the cluster. In our future 
simulations we will consider alternatives to the uniform dis- 
tribution used in the present runs. However, based on our 
current simulations it is possible to estimate the significance 
of centrally concentrating the binary population. The clus- 
ters with 50 per cent binary fractions contain approximately 
the same number of binaries within R c as a cluster with a 
10 per cent overall binary fraction in which all the binaries 
are located within R c . Given that none of our simulations 
exhibited R c growth of a magnitude sufficient to match the 
observations, a core binary fraction in excess of 50 per cent 
would be required to produce greater R c growth. The ob- 
served similarity of the d egree of mass segregat ion in the 
intermediate-age clusters Jde Griis et al.ll2002bi) render it 
unlikely that the radial distribution of binaries could vary 
so greatly between clusters. 

As in the case of the simulations in Section[3] it is natu- 
ral to ask whether our conclusions are also valid for clusters 
containing larger numbers of stars than those in our current 
simulations. The simulation of large- N clusters containing 
significant numbers of primordial binaries has only recently 
become feasible and we are currently investigating the im- 
pact of binaries on the evolution of such clusters; the results 
of this study will be presented elsewhere. The present simu- 
lations show that the disruption time for a small- TV cluster 
in a time- varying external tidal field is significantly reduced 
if it contains a large fraction of primordial binaries. It is not 
trivial to extrapolate these results to larger- TV clusters. How- 
ever, binary heating of the cluster is essentially a two-body 
process and therefore proceeds on the relaxation time scale 
£ r h- Given that a larger- N cluster with the same mean phys- 
ical density has a larger t t h, the physical rate at which the 
cluster is heated by its binary population will be reduced. 
Thus, while larger-iV clusters would obviously survive for 
longer than those in our present simulations, their physical 
rate of R c expansion would be reduced. Large variations in 
primordial binary fraction are therefore unlikely to be the 
explanation of the LMC observations. 



5 CONCLUSIONS 

The core radii of star clusters in the LMC display a strong 
trend with age, namely that the older clusters exhibit a much 
greater spread in core radius than their younger counter- 
parts. In particular, the R c of the intermediate-age clusters 
NGC 1868 and NGC 1831 differ by about a factor of three 
It ha s been suggested jElson et al.ll989l:lMackev fc Gilmorel 
2003) that the trend represents a physical evolution of the 
clusters and that the clusters all formed with similar, small 
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Figure 14. Plot of central mass density versus estimated mass 
M for the LMC cluster sample. Clusters with R c larger (smaller) 
than 4 pc are show n as op en (filled) circles. Data are taken from 
iMackev fc Gilmorei bOOffl . 



core radii which were modified by subsequent dynamical evo- 
lution. The possibility that this might be caused by varia- 
tions in the cluster IMF has been ruled out based on HST 
photom etry of a sample of eight clusters Jde Griis et alJ 
2002b). In this paper we have used A-body simulations 
of star clusters initially containing up to N = 5000 stars 
to investigate some possible alternative explanations. We 
have considered two possible dynamical effects which could 
produce significant R c evolution: time- varying external tidal 
fields and a varying population of primordial binaries. 

Throughout this paper we have analysed our sim- 
ul ations using the same procedures as were applied 
bv lMackev fc Gilmorei (120031) to the observed data. We em- 
phasised the importance of such an approach by comparing 
the evolution of the observational core radius to that of the 
half-mass radius and the density radi us. As has been pointed 
out by a number of authors (e.g. IPortegies Zwart et alJ 
l200ll) , in general there are no simple relations between these 
quantities. Thus observations of simulated data sets are es- 
sential to ensure that comparisons between observed and 
simulated data are made in a consistent manner. 

We considered clusters on circular and elliptical orbits 
of equal time-averaged radius about a point mass represent- 
ing the LMC. The clusters are scaled so that their limiting 
radii correspond to the pericentric tidal radius of the clus- 
ters on elliptical orbits. However, at the start of each simu- 
lation the clusters lie well within the tidal radius which cor- 
responds to their initial galactocentric position (apocentre 
for the clusters on elliptical orbits and the circular radius for 
the clusters on circular orbits). Initially, the clusters expand 
due to stellar mass loss but the competing effect of mass 
segregation means that _R C decreases slightly. At later times, 
two-body relaxation processes drive a further expansion of 
the clusters because the clusters are not tidally limited. This 



expansion is mirrored by a growth in R c . However, the rate 
of R c expansion is found to be identical irrespective of the 
orbit on which the cluster is moving. We conclude that dif- 
ferences in the cluster orbits cannot be responsible for the 
observed differences in R c among the LMC clusters. 

Given that the binary fractions of all but the youngest 
clusters are currently unconstrained by observations, we 
compared the evolution of the _R C of a cluster containing no 
primordial binaries with that of a cluster with a large initial 
hard binary fraction. The presence of significant numbers 
of binaries leads to increased growth of R c compared with 
the evolution of clusters containing only single stars. Binary 
fractions of 50 per cent and 25 per cent produce similar de- 
grees of i?c expansion, while a 10 per cent binary fraction 
results in noticeably less growth. Even in the most extreme 
case of a cluster with a binary fraction of 50 per cent on an 
elliptical orbit, however, R c increases by at most a factor 
of two. Not only is this smaller than the increase required 
to explain the LMC observations but the required spread 
in primordial binary fraction is already difficult to reconcile 
with the o bserved similarity of t he cluster stellar luminosity 
functions dde Griis et al.ll2002bl) . 

Our simulations have demonstrated that neither large 
variations in primordial binary fraction nor differences in 
cluster orbits about the LMC can account for the i? c -age 
relation seen in the LMC clusters. The origin of the spread 
in J? c among the intermediate-age clusters therefore remains 
an open issue. We now return briefly to the observed data 
on the LMC cluster system as a whole to see whether there 
are any further clues to the origin of the trend. 

Fig. 1141 presents evidence for a systematic difference in 
the central density of the large- R c clusters relative to the rest 
of the cluster population: these clusters have systematically 
lower central densities for a given cluster mass. This could 
be indicative of variations in star formation efficiency (SFE) 
between the clusters: low SFE leads to rapid cluster expan- 
sion and conse quently a drop in the central density (e.g. 
lGoodwirj|l997l) . Simulations of cluster expansio n due to 
resid ual gas expulsion from clusters with low SFE dGoodwinl 
do not appear to show sustained core expansion of the 
magnitude seen in the interme diate-age LMC c lusters. How- 
ever, the core radii quoted in iGoodwinl dl997l) are those of 
the mass distribution in the simulated clusters which, as 
was discussed in Section [3] may not accurately reflect the 
changes in the light distribution. We are currently investi- 
gating the long-term evolution of clusters which experience 
different degrees of primordial gas loss owing to SFEs in the 
range 0.25 — 0.60, required to expla in the range of c ore radii 
seen in the youngest LMC clusters jGoodwinlll99l) . We are 
also simulating the effects of primordial binaries on larger- 
N clusters. Of particular interest is the possibility that the 
heating effects of a large binary population could enable a 
cluster to retain a large initial R c produced by gas expulsion. 
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